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LANGEVIN  REPRESENTATION  OF  COULOMB  COLLISIONS 

IN  PIC  SIMULATIONS 

1.  Introduction 

The  development  of  particle-in-cell  (PIC)  simulation  codes’'^  was  originally 
motivated  by  work  on  the  physics  of  high-temperature  collisionless  plasmas.  In  much  of 
the  early  work,  measures  were  taken  to  reduce  the  influence  of  numerical  collisions  and 
render  the  simulation  as  close  to  collisionless  as  possible.  However,  in  recent  years  much 
interest  has  been  drawn  to  plasmas  where  collisions  play  a  central  role,  such  as  the  low- 
temperature  partially  ionized  plasmas  used  for  materials  processing.  In  these  plasmas, 
collisions  between  charged  particles  and  neutral  species  are  always  important;  indeed, 
the  gas  chemistry  that  is  central  to  the  processes  of  interest  is  often  driven  by  electron- 
neutral  (e-n)  collisions.  Thus,  methodologies  have  been  developed  for  modeling  electron- 
neutral  collisions  within  PIC  codes.  The  most  widely  used  approach  is  to  append  a 
Monte  Carlo  (MC)  step  at  which  probabilistic  scattering  occurs  between  a  randomly 
selected  charged  particle  and  a  neutral  particle.  Codes  of  this  type  are  usually  called 
PIC/MC  codes.^’^  A  related  technique  that  has  been  widely  used  in  neutral  gas 
aerodynamics  is  called  direct  simulation  Monte  Carlo  (DSMC).’^  In  most  cases,  electron- 
electron  (e-e)  scattering  has  not  been  included  in  this  type  of  code 

However,  in  high  density  processing  plasmas  such  as  electron  cyclotron  resonance 
(ECR),  helicon,  and  inductively  coupled  plasmas,  electron  densities  can  exceed  cm'^, 
and  e-e  collision  frequencies  can  exceed  those  of  e-n  collisions.  Even  in  cases  where 
pitch-angle  scattering  is  predominantly  due  to  e-n  collisions,  e-e  collisions  can  be  crucial 
in  determining  the  electron  energy  distribution  function  (EEDF).  In  many  types  of 
discharge,  the  energy  input  is  primarily  into  the  thermal  part  of  the  EEDF,  and  the  high- 
energy  tail  is  populated  primarily  by  energy  up-scattering  consequent  to  e-e  collisions. 
Electron-electron  collisions  always  drive  the  EEDF  toward  Maxwellian,  but  there  may  be 
a  competition  with  inelastic  e-n  collisions  which  deplete  the  high-energy  tail,  and  in  a 
bounded  plasma  with  escape  of  high-energy  electrons  to  the  walls.  The  high-energy  tail 
controls  atomic  excitation,  ionization,  and  to  some  extent  plasma  chemistry,  and  thus 
determines  many  of  the  properties  of  a  plasma  that  are  crucial  for  processing  applications. 
In  addition,  sheath  potentials  are  determined  by  the  competition  between  escape  of  high- 
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energy  electrons  to  the  wall  and  repopulation  of  the  high-energy  tail  by  collisions.  Thus, 
it  is  essential  that  both  e-e  and  e-n  collisions  be  modeled  accurately,  particularly  in  the 
case  of  high  density  discharges. 

We  have  recently  developed  a  2-D  axisymmetric  PIC/MC  model  of  an  ECR 
discharge  with  strongly  magnetized  electrons.*^  The  unique  feature  of  our  model  is  that 
both  electrons  and  ions  are  represented  by  particles,  but  the  electrostatic  field  is 
determined  from  the  requirement  of  quasineutrality,  rather  than  by  solving  Poisson’s 
equation.  Therefore,  plasma  oscillations  are  absent  from  the  model,  and  it  is  possible  to 
use  time  steps  many  orders  of  magnitude  longer  than  the  electron  plasma  period,  as  well 
as  spatial  gridding  much  coarser  than  the  Debye  length  Xd-  In  the  model,  elastic, 
inelastic,  and  ionizing  e-n  collisions  are  handled  with  a  Monte  Carlo  collision  scheme.  In 
this  paper,  we  describe  the  formulation  which  we  have  developed  to  include  electron- 
electron  and  electron-ion  collisions.  We  believe  it  is  suitable  for  use  in  a  wide  variety  of 
simulation  applications. 

In  DSMC  simulations  of  low-density  neutral  gases,  collisions  are  modeled  by 
picking  out  nearby  pairs  of  particles,  at  regular  intervals,  and  allowing  various  types  of 
collisions  to  occur  between  them.  The  determination  of  which  type  of  collision  occurs 
(e.g.,  elastic,  excitation,  ionization,  scattering  angle,  etc.)  depends  on  the  choice  of 
random  numbers,  with  probabilities  determined  by  the  relevant  cross-sections.  In 
.  PIC/MC  simulations  of  plasmas,  collisions  between  charged  and  neutral  particles  are 
treated  in  a  somewhat  similar  way,  but  usually  the  grid  is  used  as  an  intermediary  in  the 
collision  process,  i.e.  the  density  of  a  particular  neutral  species  is  laid  down  on  the  grid, 
and  then  the  probability  of  an  electron  colliding  with  that  species  is  proportional  to  the 
density.  The  most  straightforward  way  to  represent  e-e  collisions  in  a  PIC  code  would  be 
to  use  the  DSMC  procedure,  i.e.,  at  appropriate  time  intervals,  to  pick  out  a  number  of 
pairs  of  electrons  and  collide  them  with  the  statistics  appropriate  to  individual  electron- 
electron  collisions.  The  problem  with  this  approach  is  that  e-e  collisions  occur 
predominantly  at  long  range,  so  that  they  are  actually  a  succession  of  very  many  small 
angle  scatterings.  In  order  to  represent  individual  collisions  with  any  degree  of  accuracy, 
it  would  be  necessary  to  use  an  extremely  small  time  step.  In  fact,  even  at  a  given  instant 
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of  time,  an  electron  will  typically  be  scattering  off  many  other  electrons  simultaneously. 
Thus  it  is  numerically  inefficient,  and  really  inappropriate  physically,  to  treat  e-e 
scattering  as  a  sequence  of  MC  collisions.  Weng  and  Kushner'^  used  an  approach  rather 
more  in  the  spirit  of  plasma  PIC/MC,  where  electrons  collided  off  electron  density/energy 
distributions  laid  down  on  the  grid,  with  statistics  chosen  as  a  rough  approximation  to  the 
Coulomb  cross  section  for  individual  e-e  collisions.  Although  this  approach  has  some 
numerical  advantages  over  the  DSMC  approach,  it  still  suffers  from  the  requirement  of  an 
extremely  small  time  step  to  resolve  the  time  between  individual  collisions. 

An  alternative  approach  which  has  been  emphasized  in  the  analytical  development 
of  plasma  kinetic  theory  is  to  represent  Coulomb  scattering  through  a  Fokker-Planck 
equation.’’**^  In  the  context  of  a  PIC  simulation,  it  is  possible  to  construct  a  Langevin 
equation  (comprising  a  deterministic  friction  and  a  random  diffusive  scattering)  which  is 
entirely  equivalent  to  any  given  Fokker-Planck  equation.^®  The  dynamical  friction  and 
stochastic  diffusion  coefficients  for  the  Langevin  equation  are  represented  as  velocity- 
dependent  grid  quantities,  which  are  simply  added  to  the  macroscopic  electric  and 
magnetic  forces  acting  on  the  electrons.  Thus,  e-e  scattering  is  in  effect  represented  as  a 
scattering  of  a  single  electron  off  the  grid,  rather  than  as  a  pairwise  process. 

A  grid-based  Langevin  formalism  for  e-e  scattering  has  recently  been  developed 
by  Jones,  et  al.^’  In  this  paper,  a  key  simplification  is  made,  which  eases  the 
implementation  of  the  scheme.  This  is  the  use  of  a  velocity-independent  dynamic  friction 
Fd  and  scalar  velocity-space  diffusion  coefficient  D,  equivalent  to  treating  e-e  scattering 
as  an  isotropic  scattering  event  with  a  collision  frequency  v  that  is  independent  of  the 
velocity  of  the  test  electron.  Specifically, 


Fd=-vv, 

(la) 

D  =  2vTe/m, 

(lb) 

where  v  is  the  velocity  of  the  test  particle,  and  Te  is  the  electron  temperature.  These 
coefficients  Fd  and  D  are  consistent  with  each  other,  and  therefore  the  Langevin  equation 
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conserves  momentum  and  energy  (to  first  order  in  the  time  step),  and  correctly  drives  the 
EEDF  toward  a  Maxwellian  distribution  with  zero  mean  velocity  and  temperature  Te- 
Furthermore,  Jones  et  al^’  use  the  Spitzer  coefficients  Fd  and  D,  which  are  the  appropriate 
averages  over  an  assumed  Maxwellian  velocity  distribution  of  test  and  field  electrons. 
For  many  applications,  these  properties  may  be  sufficient.  However,  the  Coulomb 
scattering  process  is  in  fact  highly  anisotropic  (nearly  all  scattering  events  are  very  weak) 
and  strongly  velocity  dependent,  falling  off  as  v'^  for  superthermal  electrons.  As  a  result 
the  coefficients  Fd  and  D  are  greatly  overestimated  in  Eqs.  (1),  for  superthermal  electrons. 
The  time  scale  for  populating  the  tail  of  the  EEDF  is  very  much  understated  by  these 
equations,  and  if  there  is  competition  between  e-e  collisions  and  inelastic  e-n  collisions, 
the  use  of  Eqs.  (1)  can  give  a  very  inaccurate  picture  of  the  EEDF  in  steady  state. 

In  this  paper,  we  build  on  the  work  of  Jones,  et  al,^‘  to  construct  a  Langevin 
scattering  formalism  which  accurately  represents  the  multiple  small-angle  Coulomb 
scattering  process,  with  velocity-dependent  friction  and  diffusion  coefficients  derived 
from  the  actual  electron  distribution.  The  derivation  of  these  coefficients  is  reviewed  in 
Section  2.  Some  approximations  are  made  which  enormously  simplify  the  formulation 
and  reduce  the  size  of  the  data  sets  needed.  The  most  notable  of  these  is  the  assumption 
of  an  isotropic  electron  velocity  distribution  function  in  calculating  the  dynamical  friction 
Fd  and  the  diffusion  tensor  D.  In  Sec.  3  we  consider  the  application  of  the  formalism  to  a 
magnetized-electron  plasma  such  as  an  our  discharge  plasma.  In  order  to  reduce 
calculational  time,  data  complexity  and  statistical  fluctuations,  the  normalized  EEDF  used 
to  calculate  Fd  and  D  is  averaged  over  a  field  line  (but  the  actual  electron  density  at  each 
grid  cell  is  used.).  The  basic  Langevin  formulation  is  energy  and  momentum  conserving, 
but  some  of  the  approximations,  finite  time  steps,  and  statistical  fluctuations  can 
introduce  minor  deviations  from  conservation.  In  Secs.  2  and  4  we  show  how  to  restore 
exact  conservation  in  an  efficient  way.  In  Sec.  5  we  discuss  the  extension  of  the 
Langevin  scattering  formalism  to  electron-ion  scattering,  which  is  in  fact  very  much 
simpler  than  e-e  scattering.  In  Sec.  6,  we  show  the  results  of  several  computational 
exercises,  which  demonstrate  the  way  in  which  e-e  collisions  drive  the  electron 
distribution  first  toward  isotropy,  and  then  (in  the  absence  of  other  collisional  processes) 
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toward  a  Maxwellian.  We  also  show,  using  a  simplified  model  of  an  argon  plasma,  how 
the  competition  between  e-e  collisions  and  inelastic  e-n  collisions  determines  the  high- 
energy  tail  of  the  EEDF.  In  Sec.  7  the  results  are  summarized. 
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2.  General  Formulation 

A.  Review  of  the  Derivation  of  the  Fokker*PIanck  Equation  for  Coulomb  Scattering 

The  formulation  of  the  Langevin  equation  for  the  electrons  starts  with  the 
Boltzman  collision  integral  for  electron-electron  collisions: 


3f(v,t) 

at 


ee 


=  nj  d^vg^[f(v')f(v')-f(v)f(v)]. 


(2) 


where  the  tilde  indicates  the  field  electron  with  which  the  electron  of  interest  is  colliding, 
the  prime  refers  to  the  value  of  a  quantity  after  a  collision  and  unprimed  denotes  the  value 
before  the  collision,  g sv- v  is  the  relative  velocity,  n  is  the  electron  number  density, 
and  the  electron  distribution  functions  are  normalized  to  unity.  The  Coulomb  scattering 
cross  section  is^® 


da _ e^ _ 

dn"  mV  sin^0/2) 


(3) 


where  0  is  the  scattering  angle  in  the  center  of  mass  frame.  For  electron-electron 
scattering,  0  is  related  to  the  impact  parameter  b  by 


0  =  2tan  ’ 


mg^b 


(4) 


The  relative  velocity  after  the  scattering  is 
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g'  =  g  +  Ag, 


(5) 


and  of  course,  the  scattering  is  a  rotation  in  the  center  of  mass  frame,  so  g'  =  g.  If  the  3- 
coordinate  is  taken  parallel  to  g,  and  the  1,2  coordinates  perpendicular  to  g,  then 

Ag  =  g  { sin0  cos(j),  sin0  sin^,  -2  sin2(0/2) }  (6) 


and 


Vi'  =  Vi  +  Ag/2,  Vj'  =  Vj  -  Ag/2  (7) 

Because  of  the  sin'^(0/2)  dependence  of  the  scattering  cross  section  in  the 
Coulomb  potential,  the  collisions  are  dominated  by  multiple  small  angle  scattering.  In 
fact,  when  integrating  Eq.  (2)  over  0,  the  result  diverges  logarithmically  as  the  lower  limit 
of  integration  0m  approaches  zero.  This  divergence  is  resolved  by  assuming  that  the 
Coulomb  force  is  shielded  over  a  distance  of  order  the  Debye  length,  and  therefore  setting 
the  minimum  scattering  angle  0m  equal  to 


0^=2tan 


-1 


2e 


2  ^ 


mg^Xi) 


(8) 


The  terms  f(v')  and  f(v')  in  the  integrand  can  then  be  expanded  in  powers  of  Ag.  It  is 
not  difficult  to  show  that  the  only  terms  which  suffer  the  divergence  as  In  0m  are  the  terms 
involving  the  first  and  second  derivatives  of  f(v).  These  are  therefore  the  dominant  terms, 
and  it  is  appropriate  to  neglect  higher  order.  One  further  approximation  is  made:  in  the 
specification  of  0m,  Eq.  (8),  the  relative  velocity  g  of  the  pair  of  colliding  electrons  is 
replaced  by  the  thermal  average  electron  velocity  Vg.  Since  the  dependence  on  0m  is  very 
weak  (logarithmic),  the  results  are  insensitive  to  this  approximation,  which  greatly 
simplifies  the  formalism.  This  leads  to  the  standard  expression’^'*’ 
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at 


ee 


:D(v)f(v) 


(9) 


where 


^  ,  47nie^ ,  an 


(10a) 


D(v)  = 


47nie^  .  a^G 


m‘ 


avav  ’ 


(lOb) 


(11a) 


em=2tan' 


^  2e^  ^ 

vmve^^oy 


(11b) 


,  =  2j, 


H(v)  =  2|d’vj^ 


(12a) 


G(v)  =  Jd^vf(v)lv-v 


(12b) 


This  then  defines  the  Fokker-Planck  equation  for  electron-electron  scattering.  The 
coefficients  G(v)  and  H(v),  which  govern  the  diffusion  and  dynamic  friction,  are  scalar 
functions  of  the  vector  velocity  v.  Since  the  Fokker-Planck  equation  is  the  lowest  order 
expansion  of  the  Boltzmann  collision  integral  in  powers  of  it  retains  the  important 
characteristics  of  the  Boltzmann  collision  integral.  These  include  the  H  theorem  (i.e.  it 
drives  the  electron  distribution  function  to  a  Maxwellian),  as  well  as  conservation  of 
energy  and  momentum. 
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B.  Isotropic  Scatterer  Approximation 

In  typical  applications  to  PIC  codes,  it  would  be  completely  impractical  (in  terms 
of  numbers  of  particles,  computation  time,  and  statistical  fluctuations)  to  actually 
compute  the  coefficients  H(v)  and  G(v)  as  multiple  integrals,  and  then  perform  numerical 
differentiations.  However,  the  Fokker-Planck  equation  can  be  reduced  to  a  much  more 
tractable  form  by  assuming  in  Eqs.  (12)  that  the  distribution  function  f(v)  of  scatterer 
electrons  is  a  function  of  only  the  magnitude  of  the  velocity,  in  the  reference  frame  in 
which  the  electron  fluid  velocity  is  zero.  This  approximation  enormously  simplifies  and 
speeds  up  the  calculation.  In  collisional  systems,  it  is  normally  a  reasonably  good 
approximation,  since  all  electron  collisions  tend  to  isotropize  the  electron  distribution 
function.  The  thermal  part  of  the  electron  distribution  isotropizes  particularly  rapidly, 
and  e-e  scattering  of  any  electron  (even  a  fast  one)  is  normally  dominated  by  scattering 
off  thermal  electrons.  We  emphasize  that  in  the  Fokker-Planck  equation  (9),  it  is  not 
necessary  to  assume  that  the  test  particle  distribution  f(v)  is  isotropic.  It  is  reasonably 
accurate  to  follow  the  evolution  of  an  anisotropic  test  particle  distribution  while 
assuming  that  the  scattering  is  off  an  isotropic  distribution  of  field  particles.  However, 
this  approximation  does  interfere  with  the  exact  conservation  of  momentum  and  energy 
which  is  a  property  of  Eq.  (9).  In  Sec.  4  we  discuss  methods  to  insure  conservation. 

When  the  assumption  of  isotropic  scatterers  is  made,  the  integrals  over  the  polar 
and  azimuthal  angles  in  velocity  space  can  be  done  in  closed  form,  and  Eqs.  (12)  reduce 
to 


H(  v)  =  ^  dv  V  ^f  ( V )  +  8jm  J”  dv  vf  ( v ) , 


(13a) 


^f(v,.rdvv(v^ 


-H3v2)f(v) 


(13b) 
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The  velocity  derivatives  in  Eqs.  (10)  can  then  be  calculated  analytically  from  Eqs.  (13), 
which  greatly  reduces  noise  in  the  simulation.  We  find  Fd(v)  =  Fd(v)(v/v),  with 


Fd(v)  = 


47me^  ,  dH 

— 2-^1“ 
dv 


327t^ne'‘ 

2  2 
m  V 


X£dvv^f(v). 


(14a) 


Since  G  is  a  scalar  function  of  the  scalar  variable  v,  the  tensor  9^G/3v3v  is  diagonal  in  a 
coordinate  system  where  the  3-component  is  parallel  to  v.  The  only  non-zero  components 
of  the  tensor  D  are 


D33(v)  = 


47me^ 


m 


d^G 

dvridv 


47tne^  .  d^G 


327t^ne^ 

3m^ 


3cyv3  m 


dv' 


^  -^JJdvv'^f(v)  +  J“dvvf(v) 


(14b) 


D„(v)  =  D22(v)  = 


47tne^ 3^G 


47me^  ,  1  dG 
m  V  dv 


IbTt^ne'* 

3m^ 


r  dvjdvj 

X  ^JJdvv2(3v2-v2)f(v)  +  2j;dvvf(v)j, 


(14c) 


Another  type  of  approximation  is  usually  necessary,  due  to  the  fact  that  the 
number  of  particles  in  a  simulation  will  normally  be  too  small  to  calculate  the  integrals  in 
Eqs.  (14)  at  every  grid  point,  and  the  time  involved  would  be  inordinate.  Therefore,  it  is 
necessary  to  perform  some  type  of  spatial  or  temporal  averaging  in  this  step.  In  Sec.  3, 
we  show  exactly  how  we  choose  to  do  this  in  our  magnetized-electron  simulation. 

C.  Formulation  of  the  Langevin  Equation 
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In  order  to  utilize  this  formulation  in  a  PIC  code,  it  is  necessary  to  go  from  the 
Fokker-Planck  equation  to  the  Langevin  equation.  To  first  order  accuracy  in  At,  the 
Langevin  equation  in  the  form 

Av  =  FdAt  +  Q,  (15) 


is  equivalent  to  the  Fokker-Planck  equation  (9).  Here  Av  is  the  change  in  a  particle’s 
velocity,  due  to  e-e  scattering,  during  a  finite  time  step  At,  Fd  is  the  dynamical  friction, 
and  Q  is  a  random  velocity  vector  chosen  from  the  distribution 


4>(Q)  = 


_ 1 _ 

(27cAt)^^2DjiD33*'2 


f 

exp  ■ 

V 


2D„At  ; 


(16) 


However,  using  Eqs.  (14)  -  (16)  and  taking  averages  over  the  stochastic  variable  Q,  one 
can  easily  show  that  there  is  an  error  of  order  (At)^,  always  positive,  in  the  total  electron 
energy  e  after  the  collision  step,  i.e. 


<£'>-£  =  27tm(At)^J“  dvv^Fd^f(v).  (17) 

Exact  energy  conservation  (as  an  ensemble  average  over  the  stochastic  variable  Q)  can  be 
restored  by  simply  adding  a  correction  6F(v)  to  Fd(v),  specified  by  the  equation 


1  +  - 


V  > 


At  „  n  At  _  , 


(18) 


For  most  purposes,  it  is  more  than  adequate  to  insure  that  the  Langevin  equation  • 

conserves  energy  to  second  order  in  At,  which  will  hold  if  5F(v)  is  given  by  the  much 

simpler  approximate  form  of  Eq.  ( 1 8),  * 
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(19) 


6F(v)  =  -^Fd2(v). 

We  note  in  passing  a  surprising  property  of  Eq.  (14a):  for  a  test  electron  with 
speed  V,  scattering  off  electrons  with  speed  v  >  v  does  not  contribute  to  the  friction 
coefficient  Fd(v),  This  is  a  peculiarity  of  Coulomb  scattering  off  an  isotropic  distribution 
of  scatterers.  As  a  result,  the  friction  is  very  small  for  electrons  with  velocity  much  less 
than  the  thermal  velocity  Ve,  and  if  the  EEDF  happens  to  be  hollow,  with  v  >  Vmin  for  all 
electrons,  then  an  electron  with  v  =  Vmin  will  feel  no  friction  at  all.  At  first  sight,  it  seems 
paradoxical  that  such  electrons  will  not  be  driven  toward  a  distribution  centered  about  the 
mean  electron  velocity  Ue.  However,  what  actually  happens  is  that  a  low  velocity  electron 
will,  on  the  average,  diffuse  up  toward  speed  Ve  (the  diffusion  coefficients  have 
contributions  from  scattering  off  electrons  with  v  >  v  ,  and  do  not  go  to  zero  for  small  v), 
and  then  feel  a  strong  friction  that  tends  on  the  average  to  center  it  about  14. 

3.  Application  to  a  Magnetized  Plasma 

A.  Langevin  Scattering  Formulation  for  Magnetized  Electrons 

This  electron-electron  scattering  model  was  developed  primarily  for  use  in  our 
recently  developed  2D-3v  (axisymmetric)  simulation  model  of  ECR  processing 
plasmas.’^  In  this  code,  both  the  electrons  and  ions  are  represented  as  simulation 
particles,  subject  to  a  strong  external  magnetic  field,  self-consistently  determined 
electrostatic  fields,  and  collisions.  The  ions  are  not  typically  strongly  magnetized;  hence 
their  trajectories  are  followed  in  full  detail  in  two  spatial  and  three  velocity  coordinates. 
However,  the  electrons  are  strongly  magnetized,  and  can  be  regarded  as  firmly  attached  to 
a  given  magnetic  field  line.  [Because  of  the  axisymmetry,  all  drifts  are  azimuthal,  and 
thus  do  not  affect  the  electron  trajectory  in  the  r-z  plane  of  the  simulation.]  Thus,  we  use 
the  actual  field  lines  as  one  set  of  elements  for  a  curvilinear  grid,  allowing  us  to  specify 
the  position  of  an  electron  by  the  field  line  number  to  which  it  is  permanently  attached, 
and  a  single  axial  coordinate  z  giving  its  location  on  the  field  line.  In  velocity  space,  the 


11 


code  follows  the  parallel  velocity  vn  of  each  electron,  and  the  magnitude  of  the 
perpendicular  velocity  v^,  but  there  is  no  need  to  follow  the  phase  of  the  perpendicular 
motion.  In  practice,  we  follow  the  magnetic  moment  p  =  mvx^/21BI,  which  is  constant 
during  the  interval  between  collisions,  rather  than  propagating  vj.  itself.  The  unique 
feature  of  our  simulation  is  that  Poisson's  equation  is  not  used,  but  rather  the  electric  field 
is  obtained  from  the  quasi-neutrality,  much  like  the  procedure  that  is  used  in  fluid 
simulations.  This  allows  the  simulation  to  avoid  inverse  plasma  frequency  time  scales 
and  Debye  length  scales.  For  both  electrons  and  ions,  the  relevant  length  scale  is  the 
macroscopic  length  scale,  and  the  minimum  time  scale  is  this  length  scale  divided  by  the 
electron  thermal  velocity.  This  simulation  scheme  speeds  up  the  calculation  by  several 
orders  of  magnitude  (for  the  usual  ECR  reactor)  as  compared  to  a  standard  particle  in  cell 
code  with  the  electrostatic  field  calculated  from  Poisson's  equation.  The  details  of  the 
simulation  scheme  and  some  results  are  given  elsewhere.’^ 

In  applying  the  e-e  scattering  formalism  to  this  type  of  guiding  center  electron 
model,  it  is  necessary  only  to  specify  what  vu'  and  \x  are  after  the  collision,  given  vh  and 
vi  before  the  collision.  To  transform  the  Langevin  scattering  results  to  these  variables,  it 
is  convenient  to  use  a  coordinate  system  with  the  3-coordinate  parallel  to  v,  the  electron’s 
velocity  before  collision,  the  1 -coordinate  normal  to  v  but  in  the  B-v  plane,  and  the  2- 
coordinate  normal  to  both  B  and  v,  as  shown  in  Fig.  1.  Let  a  be  the  angle  between  B  and 
V,  so  that  tan  a= vx/ V||.  After  scattering,  the  new  velocity  v'  is  given  by 


v,'  =  Qi, 

(20a) 

V2'  =  Q2. 

(20b) 

V3'  =  V  +  FjAt  +  Q3. 

(20c) 

Transforming  back  to  vn'  and  vx',  we  find  from  Fig.  1  that 

V||'  =  va'cos  a  -  v/  sin  a  =  ( V  +  FjAt  +  Q3)cos  a  -  Qi  sin  a,  (21a) 

vx'^  =  (vs'  sin  a  +  v/  cos  a)^  +  \2^  =  [(v  +  F^At  +  Q3)  sin  a  +  Qi  cos  a]^  +  Q2^’  (21b) 
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Equations  (21),  together  with  Eqs.  (14)-(16),  give  the  basic  Langevin  formulation  for  the 
e-e  scattering  in  a  magnetized  system. 

B.  Data  Structure  and  Numerical  Considerations 

By  assuming  that  the  distribution  of  scatterers  is  isotropic,  we  have  reduced  the 
scattering  coefficients  Fd(v),)  Dii(v)  and  D33(v)  to  grid-dependent  quantities  that  also 
depend  on  the  magnitude  of  the  electron  velocity.  However,  statistical  fluctuations 
incident  to  the  finite  number  of  simulation  particles  would  make  it  virtually  impossible  to 
actually  compute  these  velocity-dependent  coefficients  by  performing  the  velocity 
integrals  of  Eq.  (14)  at  every  grid  point.  (It  would  also  be  inordinately  time-consuming.) 
For  example,  in  a  two  dimensional  PIC  simulation  with  a  100x100  grid,  even  with  a 
million  particles  there  are  only  100  particles  per  cell,  leading  to  fluctuations  at  least  on 
the  order  of  10%  (and  even  worse  for  v-dependent  quantities).  Clearly,  some  spatial 
and/or  temporal  averaging  is  necessary.  In  our  magnetized-electron  code,  there  is  an 
obvious  way  to  do  this.  The  normalized  velocity  distribution  fj(v)  is  constructed  for  all  of 
the  electrons  on  field  line  j,  and  is  then  used  to  calculate  the  velocity  integrals  in  Eqs. 
(14).  However,  the  density  n  used  in  Eq.  (14)  is  taken  to  be  the  local  electron  density  on 
the  grid.  Since  electrons,  in  their  orbits,  rapidly  sample  an  entire  field  line,  this  should  be 
a  very  good  approximation.  It  is  also  very  efficient  numerically,  since  the  integral 
quantities  in  (14)  can  simply  be  accumulated  at  the  same  time  that  the  particle  densities 
are  laid  down  on  the  grid. 

In  the  next  section,  we  show  how  to  build  exact  energy  and  momentum 
conservation  into  the  scattering  formalism.  If  the  conservative  form  (22)  or  (23)  is  used, 
the  only  limits  imposed  on  the  scattering  time  step  At  are  those  necessary  to  insure 
accuracy.  Thus,  At  should  be  no  more  than  a  fraction  of  the  e-e  collisional  relaxation 
time.  In  many  cases,  there  will  be  stronger  constraints  imposed  by  other  aspects  of  the 
simulation.  For  example,  if  e-e  scattering  is  competing  with  inelastic  electron-neutral 
scattering,  then  At  should  be  no  larger  than  the  characteristic  time  for  the  latter  process. 
In  a  bounded  plasma,  sheath  potentials  may  be  determined  by  competition  between 
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escape  of  high-energy  electron  through  the  sheath,  and  replenishment  of  high-energy 
electrons  via  e-e  collisions.  Then  At  must  be  no  larger  than  the  time  step  used  for 
electron  escape.  In  other  cases,  the  limit  on  At  may  arise  from  the  characteristic  time  for 
electrons  to  transport  spatially  from  one  region  to  another. 


4.  Energy  and  Momentum  Conservation 

The  Fokker-Planck  equation  (9),  with  coefficients  from  Eqs.  (10)  -  (12),  exactly 
conserves  momentum  and  energy.  However,  the  numerical  implementation  of  e-e 
scattering,  as  described  in  Secs.  2  and  3,  may  suffer  small  deviations  from  momentum 
and  energy  conservation.  Errors  in  energy  conservation  are  particularly  troublesome, 
since  Lemons,  et  al^^  have  shown  that  they  can  lead  to  systematic  (non-random)  drifts  in 
total  energy  that  can  become  substantial  over  long  times.  These  types  of  effects  could 
lead  to  significant  errors  in  EEDF,  and  therefore  in  ionization  fraction,  chemical  make-up, 
and  other  such  properties  of  the  plasma.  Momentum  conservation  errors  may  also  be  of 
concern,  if  they  interfere  with  the  calculation  of  electric  currents  to  sufficient  accuracy. 

Momentum  and  energy  non-conservation  can  occur  for  two  reasons.  First,  if  the 
electron  velocity  distribution  f(v)  is  not  exactly  isotropic,  the  assumption  of  an  isotropic 
distribution  f(v)  for  the  scatterers  is  inconsistent  with  the  conservation  laws,  even  if  the 
mean  electron  velocity  is  zero.  Obviously,  the  situation  will  be  worse  if  Ue  is  non-zero 
and  the  isotropic  distribution  f(v)  is  prepared  in  a  frame  of  reference  other  than  that  of  u*. 
In  Sec.  3,  we  have  chosen  to  calculate  f(v)  in  the  lab  frame  for  all  the  electrons  on  a  given 
field  line,  so  if  Ue  is  non-uniform  along  the  field  line,  the  scattering  formalism  clearly  will 
not  conserve  electron  momentum  locally.  Secondly,  the  diffusion  part  of  the  Langevin 
equation  involves  the  choice  of  random  velocity  increments  Q.  On  the  average,  these 
increments  will  conserve  energy  and  momentum,  but  given  the  finite  number  of 

simulation  particles  N  at  any  grid  cell,  errors  of  the  order  of  Vn  can  always  occur. 

Fortunately,  these  non-conservation  effects  are  small,  particularly  if  Ug «  Ve,  as  is 
the  case  in  our  ECR  plasma  simulations  and  many  other  typical  situations.  It  is  then 
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particularly  easy  to  make  corrections  that  adequately  restore  the  conservation  laws.  In 
general,  one  can  restore  momentum  and  energy  conservation  locally  by  using  a 
renormalization  procedure  described  by  Lemons,  et  al.^^  We  take  note  of  the  electron 
fluid  velocity  Ucii(z)  and  the  electron  temperature  Te(z)  on  field  line  j  before  the  e-e 
scattering  step,  and  the  values  Ueii'(z)  and  Te'(z)  after  scattering.  In  general,  they  will  be 
slightly  different.  We  then  reset  the  velocity  v/  of  electron  n  located  at  z,  according  to 
the  formula 


In  practice,  we  find  that  Eq.  (22)  is  overkill.  Small  random  errors  in  electron 
momentum  conservation  during  e-e  collisions  are  usually  unimportant,  since  collisions  of 
electrons  with  neutrals  and/or  control  the  electron  current.  Furthermore,  it  is  usually 
sufficient  to  insure  energy  conservation  globally  over  some  large  area,  in  our  case  over  a 
field  line,  since  the  velocity  integrals  in  Eq.  (14)  are  performed  as  an  average  over  a  field 
line.  Thus  we  use  the  very  simple  velocity  renormalization 


(23) 


where  Wj  is  the  total  kinetic  energy  of  all  the  electrons  on  field  line  j  before  the  e-e 
collision  step,  and  Wj'  is  the  same  quantity  after  the  collision  step. 


5.  Electron-Ion  Scattering 

The  Fokker-Planck  equation  for  electron-ion  (e-i)  scattering  is  derived  in  exactly 
the  same  way  as  Eqs.  (9)  -  (12).  The  only  difference  in  the  results  is  that  Eqs.  (12)  are 
replaced  with 
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(24a) 


nij  J  |v-v| 


G(v)  =  Zi2Jd^vfi(v)|v-v|, 


(24b) 


where  Zje  is  the  ion  charge  and  fi(v)  is  the  ion  velocity  distribution. 

The  rate  of  energy  exchange  between  electrons  and  ions  is  down  by  order  me/nii, 
which  makes  it  negligible  for  many  purposes.  If  we  neglect  energy  exchange,  and  treat  e- 
i  scattering  as  essentially  just  pitch-angle  scattering  of  the  electrons  off  infinitely  massive 
ions,  then  the  formalism  becomes  particularly  simple.  It  is  then  appropriate  to 
approximate  the  ion  velocities,  which  are  always  small  compared  to  Ve,  as  zero,  so  that 
Eqs.  (24)  reduce  simply  to 

H  =  ^,  (25b) 

V 

G  =  Zi\.  (25b) 

According  to  Eqs.  (14),  the  dynamical  friction  coefficient  is 


Fd=- — 
m  V 


(26a) 


and  the  diffusion  coefficients  are 


D33(v)  =  0, 


(26b) 


Di,(v)  =  D22(v)  = 


4nne%^ , 

2  ^ 
m  V 


(26c) 
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For  scattering  of  electrons  off  infinitely  massive  ions,  momentum  conservation  is 
not  a  consideration.  Electron  energy  should  be  conserved  exactly  in  every  collision,  so 
the  simplest  procedure  is  to  not  use  the  dynamical  friction  from  Eq.  (26a),  but  simply  to 
specify  Avs  so  as  to  insure  exact  energy  conservation: 

(v+Av3)V  Qi  V  Qa  =  v^  (27) 

where  Qi  and  Q2  are  the  stochastic  increments  to  the  velocity  components  normal  to  v, 
chosen  from  the  distribution  (16).  If  we  neglect  second  order  in  Qi/v  and  Q2/V,  Eq.  (27) 
becomes  simply 


AV3 


Qi^  +  Q2^ 

2v 


(28) 


and  in  the  case  of  magnetized  electrons,  the  electron  velocity  components  parallel  and 
perpendicular  to  B,  from  Eqs.  (21),  become 


/ 

V- 

V 


2v  J 


cosa-QjSina, 


(29a) 


V- 


2v 


^2 


sina  +  QjCosa 


+  Q2 


(29b) 


If  it  is  important  to  calculate  energy  transfer  between  the  electrons  and  ions,  it  is 
easy  to  modify  the  formalism  to  include  this  effect. 

6.  Computational  Examples 
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In  this  section,  we  present  some  simple  computational  examples  as  test  cases  for 
our  formulation  of  e-e  scattering.  We  consider  a  ID  system  with  periodic  boundary 
conditions.  (In  the  context  of  our  2D  ECR  plasma  code,  this  could  be  thought  of  as  a 
single  magnetic  field  line,  with  uniform  magnetic  field.  However,  the  calculations  in  this 
section  do  not  include  the  end  losses  that  would  occur  in  a  bounded  plasma.)  All  initial 
conditions  (e.g.,  density,  distribution  function)  are  spatially  uniform.  The  system  length  is 
35  cm  and  the  cell  size  is  Az  =  1  cm.  We  use  15,000  macroparticles  to  represent  the 
electrons.  The  electrons  are  scattered  according  to  Eqs.  (15)  -  (18)  at  intervals  At  =  8x10'^ 
sec. 

A.  Approach  to  Equilibrium 

In  this  case,  we  consider  the  evolution  of  the  electron  distribution  from  a  very 
anisotropic  and  non-Maxwellian  initial  condition,  with  e-e  scattering  the  only  physical 
process  represented  in  the  simulation.  The  initial  distribution  represents  two  cold  counter 
streaming  electron  beams. 


f(vii,  vx,  t=0)  =  j  [6(V||  -  Vo)  +  6(vii  +  Vo)]  5(vx),  (30) 

with  beam  energy  jmvo^  =  4  eV  The  plasma  density  is  cm'^.  Since  every  cell  is 

identical,  this  is  simply  a  point  problem  from  a  fundamental  point  of  view.  However,  it  is 
still  useful  to  think  in  terms  of  a  simulation  of  all  the  points  along  a  field  line,  since  this  is 
the  way  the  data  structures  and  statistical  properties  of  the  simulation  are  organized. 

We  recall  that  the  diffusion  coefficients  of  Eqs.  (14)  decrease  rapidly  with  particle 
speed  V,  so  that  one  expects  the  approach  to  equilibrium  to  proceed  rapidly  for  electrons 
in  the  low-energy  (thermal)  range,  and  more  slowly  in  the  high-energy  tail.  The  standard 
estimate’®  for  the  relaxation  time  is  trei  =  6xl0  ’  sec.  Fig.  2  shows  plots  of  the  reduced 
electron  distribution  functions  f(vii)  and  f(Vi) ,  at  four  different  times.  In  these  plots,  the 
abscissa  is  chosen  to  be  £«  s  imeVn^  or  ei  s  jmeVi^.  In  Fig.  2a,  at  the  early  time  t  = 

5x10  ®,  the  distribution  functions  still  show  the  cold  two-beam  structure.  Figure  2b 
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shows  the  distribution  functions  at  t  =  IxlO'^sec.  By  this  time,  f(V||)  has  become  single- 
peaked  and  is  fairly  close  to  Maxwellian  in  the  energy  range  up  to  6  eV,  but  at  higher 
energies  the  distribution  falls  off,  and  anisotropy  is  increasingly  evident.  Figure  2c  shows 
the  distributions  at  t  =  2x10'^  sec.  By  this  time,  both  f(vii)  and  f(vx)  are  close  to 
Maxwellian,  but  still  at  slightly  different  temperatures.  Fig.  2d  shows  that  at  t  =  6x10’’ 
sec  (the  predicted  relaxation  time),  the  distribution  functions  are  isotropic  and 
Maxwellian  over  their  entire  energy  range. 

B.  Balance  Between  Heating,  e*e  Collisions,  and  Inelastic  Collisions 

In  this  example  we  model,  in  a  very  simplified  way,  the  combined  effect  of 
several  processes  that  occur  in  an  ECR  discharge:  plasma  heating,  electron-electron 
collisions,  and  electron  energy  loss  due  to  ionizing  collisions.  In  the  model,  electrons 
with  energy  e  <  3  eV  are  heated  every  time  they  pass  a  “resonant  zone”.  The  heating  is 
implemented  by  giving  each  electron  a  velocity  kick  each  time  it  passes  by  the  position  z 
=  3  cm,  with  the  velocity  increment  Av  chosen  randomly  from  a  Gaussian  distribution 
with  mean  value  |m(Av)^  =  1  eV.  We  also  include  electron  energy  loss  due  to  “ionizing 

collisions”  with  neutral  atoms.  (However,  we  do  not  create  new  electrons  when  one  of 
these  “ionizing  collisions”  occurs.  Since  we  are  simulating  a  closed  system  with  no 
particle  losses,  this  would  preclude  the  attainment  of  steady  state.)  The  ionization  cross 
section  for  Ar  is  used,  as  given  by  Tachibana^^.  This  cross  section  increases  from  about 
10'’^  cm^  just  above  the  ionization  threshold  =  15.76  eV,  to  a  maximum  of  3.9xl0'® 
cm^  at  60  eV.  The  neutral  gas  pressure  is  taken  to  be  5  mTorr.  Each  electron  loses 
exactly  15.76  eV  of  energy  when  it  is  scattered.  Electron-impact  excitation  collisions  are 
omitted  from  this  simple  model,  even  though  they  do  represent  a  significant  energy  loss 
mechanism  in  a  real  gas. 

Figure  3a  shows  the  equilibrium  electron  energy  distribution  function  for  a  case 
with  plasma  density  De  =  10  cm* .  At  this  value  of  ne,  electron-electron  scattering  is 
weak,  and  ionization  energy  losses  deplete  the  tail  of  the  distribution  function  for  energies 
above  Ciz.  Figure  3b  shows  the  EEDF  for  a  plasma  of  density  lO”  cm*^.  Here,  e-e 
scattering  is  strong  enough  to  drive  the  electron  distribution  to  Maxwellian  in  the  regime 
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below  Eiz,  and  to  significantly  replenish  the  distribution  above  the  ionization  threshold. 
Figure  3c  shows  the  EEDF  for  plasma  density  10*^  cm'^.  In  this  case,  e-e  scattering  is 
easily  strong  enough  to  redistribute  energy  from  the  heating  region  e  <  5  eV  to  the  tail 
region,  and  the  equilibrium  distribution  is  very  nearly  Maxwellian  over  the  entire  energy 
range. 


7.  Conclusions 

The  Langevin  equation  can  be  used  to  formulate  electron-electron  and  electron- 
ion  collisions  in  a  probabilistic  manner  analogous  to  the  Monte  Carlo  treatment  that  is 
often  applied  to  electron-neutral  or  neutral-neutral  scattering.  The  difference  is  that  e-e 
and  e-i  collisions  are  very  frequent  and  very  weak,  so  that  the  Langevin  equation 
represents  the  net  probabilistic  effect  of  very  many  small-angle  scatterings.  Therefore, 
time  steps  can  be  long  compared  to  the  time  scale  for  interactions  between  particular  pairs 
of  charged  particles.  The  exact  form  of  the  Langevin  equation  is  well  known,  but  is 
impractical  for  numerical  applications,  due  to  the  need  for  very  large  numbers  of 
simulation  particles,  extensive  data  structures,  and  burdensome  computations.  We  have 
used  several  simple  and  well-justified  approximations  to  reduce  the  formulation  to  a 
manageable  and  efficient  form.  We  have  also  provided  simple  procedures  for  insuring 
that  momentum  and  energy  are  conserved  in  the  numerical  implementation.  The  general 
formulation  is  applicable  to  either  unmagnetized  or  magnetized  electrons,  and  in  the  latter 
case  we  have  expressed  the  results  specifically  in  terms  of  the  velocity  components 
parallel  and  perpendicular  to  B. 
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Fig.  1  —  Geometry  for  the  transformation  from  coordinate  3  (parallel  to  the  velocity  v  before  scattering),  1  (normal 
to  V  but  in  the  B-v  plane),  2  (normal  to  both  B  and  v)  to  coordinates  ||  (parallel  to  B)  and  1  (normal  to  B). 
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Fig.  2  (Continued)  —  Reduced  electron  distribution  functions  f(V|)  (solid  curve  and  f(Vj^)  (dashed  curve)  at  times 
(a)  5  X  10"®  sec,  (b)  1  x  10"’  sec,  (c)  2  x  10”’  sec,  and  (d)  6  x  10”’  sec. 


0  5  10  15  20 

Energy  (eV) 

(a) 


Fig.  3  —  Electron  energy  distribution  functions  for  Ar  at  pressure  5  mTorr,  after  steady  state  has  been  reached. 
Included  are  a  simple  model  of  bulk  electron  heating,  e-e  collisions,  and  electron  energy  losses  to  ionizing 
collisions.  Excitations,  and  electron  creation  and  loss  are  not  included,  (a)  n^  =  10'°  cm"\  (b)  n^  =  10"  cm■^ 
(c)  tie  =  10'^  cm"^. 
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Fig.  3  (Continued)  —  Electron  energy  distribution  functions  for  Ar  at  pressure  5  mTorr,  after  steady  state  has  been 
reached.  Included  are  a  simple  model  of  bulk  electron  heating,  e-e  collisions,  and  electron  energy  losses  to  ionizing 
collisions.  Excitations,  and  electron  creation  and  loss  are  not  included,  (a)  n^  =  10'°  cm”^.  (b)  n^  =  10"  cm”^. 
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Fig.  3  (Continued)  —  Electron  energy  distribution  functions  for  Ar  at  pressure  5  mTorr,  after  steady  state  has  been 
reached.  Included  are  a  simple  model  of  bulk  electron  heating,  e-e  collisions,  and  electron  energy  losses  to  ionizing 
collisions.  Excitations,  and  electron  creation  and  loss  are  not  included,  (a)  n^  =  10’°  cm”°.  (b)  n^  =  10"  cm”^ 
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